Recent Advances in the Numerical Simulations of 
Binary Black Holes 

Pedro Marronetti^, Wolfgang Tichy^ 

^ Physics Department, Florida Atlantic University, Boca Raton, FL 33431-0991 
T-H E-mail: pmarrone@fau.edu 

o 

^_^ Abstract. Since the breakthrough papers from 2005/2006, the field of numerical 
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^^ relativity from the category of "really-hard-problems" to the realm of "things-we- 

^^ know-how-to-do" . Simulations of binary black holes in circular orbits, the holy grail 

of numerical relativity, are now tractable problems that lead to some of the most 

'~~^ spectacular results in general relativity in recent years. We cover here some of the 

Q-( latest achievements and highlight the field's next challenges. 
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^ 1. Introduction 

^ Whether black holes (BH), defined as a family of solutions to the vacuum Einstein 

'k>( field equations that present a one-way membrane (i.e., the event horizon) that causally 

H separates two regions of spacetime, exist or not was a matter of speculation for the 

best part of the 20th century. In the past decade, however, astronomical observations 
placed them as the most promising models for objects detected in X-ray binaries (with 
sizes of a few to tens of solar masses) and for the supermassive entities at the center of 
most galaxies (with up to billions of solar masses) [1]. While they can be fully specified 
by their mass, angular momentum and charge, only the former two are of astrophysical 
relevance since charged BHs are quickly neutralized by free charges found in their vicinity 
(i.e, from accretion disks, interstellar plasma, etc.). 

While individual BHs are extremely interesting objects on their own, when they 
couple in binary black holes (BBH) they become more so by forming one of the 
most promising sources for gravitational wave detectors [2]. The two-body problem 
in general relativity is unsolved in the sense that in this framework we do not have 
solutions analogous to the Keplerian curves of Newtonian gravity. Due to the emission of 
gravitational waves, the binary loses energy and angular momentum tightening the orbit 
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until finally the two holes merge to form a single larger distorted BH. This still emits 
gravitational waves until it reaches the state of a stationary rotating BH. Until recently, 
the most detailed predictions existed only in two disjoint regimes where approximation 
schemes could be used. As long as the two objects are far apart, post-Newtonian 
calculations [3] can give highly accurate approximations for the orbital motion and the 
gravitational waves emitted by the binary. When the two BHs get closer, the post- 
Newtonian expansion becomes more and more inaccurate and eventually breaks down. 
Later in the evolution of the binary when the two BHs have already merged, the close 
limit approximation [4, 5, 6] can be used to model the ring-down of the remaining 
single BH. However, in the intermediate region that covers the last couple of orbits 
and up until the point when a single BH forms, the full non-linear Einstein equations 
have to be solved and computer simulations are used to obtain numerical answers [7]. 
Such simulations are a hard and challenging problem. A problem that with the first 
generation of interferometric gravitational wave detectors fully operational and on their 
way to their second generation, has become increasingly important. As sensitive as 
these instruments are, it seems likely that experimental detection will hinge on detailed 
and accurate theoretical predictions. Such predictions would help distinguish the weak 
signature of real astrophysics among the various types of noise known to be present in 
the apparatus. 

Enter numerical relativity. The field of numerical relativity is relatively young since 
its development has been closely paired to the growth of the power and availability of 
computers. While some pioneer work was done as early as the mid 1960s, it wasn't 
until the 1970s that the first results of BBH spacetime evolution were obtained [8, 9]. 
These corresponded to axi-symmetric BH collisions. Researchers soon found that the 
situation gets much more difficult when dealing with problems without a high degree of 
symmetry as in the case of BBH in semi-circular orbits. To say that these simulations 
are a hard and challenging problem is a gross understatement. Until recently they 
tended to fail after a very short time due to instabilities which resulted in exponentially 
growing run-away solutions. Fortunately, since the breakthrough papers of Pretorius 
[10], Campanelli et al. [11] and Baker et al. [12] tremendous progress has been achieved 
and it is now possible to evolve BBHs through many orbits and the subsequent merger 
and ringdown phases. We will cover here some of the most spectacular discoveries 
regarding BBH made through the solution of the full set of Einstein field equations. 
Many of the foundations of the subjects touched here are described in the excellent 
books by M. Alcubierre [13] and T. Baumgarte and S. Shapiro [14]. Quantities in this 
article are expressed Geometrized units, where G = c= 1. 
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2. Binary Black Holes: Recent Developments and Discoveries 

2.1. Orbital Hang-up and Naked Singularities 

Relativists have speculated for almost a century about the possibility of finding naked 
singularities (i.e., those that do not hide behind event horizons). Penrose in his famous 
Cosmic Censorship conjecture stated that no naked singularities exist apart from the 
one associated with the Big Bang. Choptuik showed in 1997 [15] that under certain 
very particular set of conditions "non-generic" naked singularities could be created |. 
Thus, the search for mechanisms (natural or artificial) that would lead to the formation 
of naked singularities continues. Different astrophysical scenarios such as rotating BHs 
accreting mass-energy, hypermassive neutron star collapse, hypernovae, high-speed BH 
collisions, etc., have been studied with the goal of forming stationary compact objects 
with a normalized angular momentum higher than critical [S/M"^ > 1) which would 
indicate the presence of a naked singularity. So far all these studies have lead to negative 
answers. One particular scenario is that of BBH mergers whose total normalized angular 
momentum at some point during the inspiral is larger than the critical value. Campanelli 
et al. [16] studied first this situation by considering two BBH systems composed by 
two equal-mass BHs positioned at the same separation and with dimensionless spin 
of magnitude s/m"^ = 0.757. However, in one case the spins were aligned with the 
orbital angular momentum (leading to an initial angular momentum of Ji/M^ = 1.18) 
and in the other they were anti-aligned (now with Ji/Mf = 0.57). Figure 1 shows these 
configurations. Their study showed that the dynamics of the final orbits are significantly 
affected by the orientation of the spins. The aligned case required three orbits before the 
merger, radiating ~ 7% of the total energy, while the anti-aligned case radiated only 
~ 2%, merging in less than one orbit. In both cases the spin of the final BH was less 
than critical, indicating that forming naked singularities by the way of merging BHs may 
not be possible. This conclusion has been confirmed so far by several groups that have 
studied mergers with different masses and spins of arbitrary magnitude and direction. In 
particular, the formulas that predict the final mass and spin of a BBH merger described 
in Section 2.4, albeit approximate and only accurate for BHs with comparable mass, 
fail to predict a region of parameter space for which BHs with higher-than-critical spin 
could be formed. 

2.2. Black Hole Kicks and Anti-kicks 

With the exception of highly symmetric situations such as head-on collisions of equal 
mass non-spinning BHs, generic BBH will radiate linear and angular momentum. This 
effect has important consequences for astrophysics. It means that any BBH with unequal 
masses will acquire a kick due to gravitational radiation during the inspiral and merger 
of the system. The magnitude of the resulting recoil is important in a variety of 

I Simply put, "non-generic" singularities are such that arbitrarily small perturbations would send them 
back inside event horizons. 
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Figure 1. At the top (right) we have a BBH system with components whose spins 
are anti-aligned with the orbital angular momentum, resulting in a total angular 
momentum less than critical. At the bottom (right) the same BHs are flipped to align 
the spins with the orbital angular momentum, resulting in a total angular momentum 
higher that critical. Numerical simulations show that the lower system will take a 
longer time to merge, emitting the excess of angular momentum that would have 
otherwise produce a naked singularity. The numerical values are those of Campanelli 
et al. [16]. 



z/M 



x/M 




Figure 2. Trajectories of the black holes in a BBH simulation. In this case, the 
individual BHs have arbitrary spins which force the orbital plane to precess and the 
final BH to acquire a recoil velocity ("kick") represented by the straight line anti- 
parallel to the z-axis. 
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astrophysical scenarios, such as the cosmological evolution of supermassive BHs or the 
growth and retention of intermediate-mass BHs in dense stellar clusters. For a binary in 
almost circular orbit, the direction of the instantaneous linear momentum flux rotates 
in the orbital plane with the angular velocity of the system. Thus, when the binary goes 
through one orbital period, the average linear momentum flux will be close to zero. The 
only net effect comes from the fact that the inspiral orbits are not perfect circles. Most 
of the kick is accumulated during the last orbit and subsequent plunge of the two holes, 
when the motion is no longer quasi-circular and the averaged linear momentum flux is 
much larger. Several analytical estimates of the kick velocity have been published in 
recent years [17, 18, 19, 20]. All these estimates have been derived using approximations 
(such as post-Newtonian theory) which break down during the last orbit missing the 
strongest contribution of the gravitational radiation to the kick velocity. On the other 
hand great progress has been made in the past couple of years with full GR calculations 
(see [14] and references therein). It is now possible to calculate the kick velocity for any 
initial configuration using numerical simulations. 

For non-spinning BHs the maximum kick velocity occurs for a mass ratio of 
q = mi/m2 = 2.77 and has a magnitude of 175 km/s [21]. For equal mass BHs 
a kick occurs only if the initial BHs are spinning. If the spins of the initial BHs 
are aligned or anti-aligned with the orbital angular momentum the maximum kick 
velocity has a magnitude of 448 km/s [22] directed in the orbital plane. This maximum 
occurs for maximally spinning BHs with spins that are anti-aligned with each other. 
However, much larger kicks are possible for equal-mass binaries with anti-aligned initial 
spins in the orbital plane. Simulations have yielded kicks as high as 2500 km/s [23]. 
Extrapolations from these simulations by Campanelli et al. [24] predict kick velocities 
as high as 4000 km/s. Kicks on the order of 1000 km/s seem to be a generic feature even 
if the initial spins are not exactly in the orbital plane and not exactly anti-aligned [25] . 
Some of the latter kicks are so high that the final BH would escape both from dwarf 
elliptical and spheroidal galaxies (with typical escapes velocities of below 300 km/s) 
and from giant elliptical galaxies (2000 km/s) [26]. The question thus arises how often 
such high kick velocities occur. After all, most galaxies that have undergone mergers 
seem to retain supermassive BHs at their centers. The most likely answer [27] is that in 
real galaxy mergers torques from accreting gas align the spins with the orbital angular 
momentum, reducing the maximum velocity to a few hundreds km/s. 

If one monitors the velocity of the center of mass during the inspiral and merger, one 
finds that the velocity continuously changes direction while slowly growing in magnitude 
during the inspiral. The final black hole directly after merger is at first distorted and still 
radiates while settling down to a more symmetrical stationary state. In many cases this 
radiation is directed such that center of mass speed decreases again. This phenomenon 
is known as the anti-kick [28, 29, 30]. 
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2.3. Spin Flips 

Black holes with masses ranging from millions to several billions of solar masses are 
colloquially known as supermassive and are found at the center of most galaxies. Of 
particular interest are those found in quasars and in galaxies with Active Galactic 
Nuclei (AGN) since they are believed to be the engine behind powerful jets that could 
reach kiloparsec scales. These jets are launched perpendicularly to the inner part of an 
accretion disk that surrounds the central BH and, barring a highly dynamical situation, 
they should be aligned with the BH axis of rotation (Bardeen-Petterson effect). Radio 
observations have shown that some galaxies have jets that at some point in their past 
changed abruptly directions. A class known as X-shaped radio galaxies [31, 32] show 
the most dramatic cases, where the radio lobes in the inner region (younger) are shifted 
by almost ninety degrees with respect to those in the outer zones (older). This sudden 
change of jet orientation is associated with a change in the direction of the spin of the 
BH. Several theoretical models have been proposed to explain these phenomena that 
involve backflow from the main lobes and instabilities of the accretion disk. However, 
the most widely accepted explanation is the model that invokes a BBH merger. The 
idea is that the original AGN galaxy merges with a neighboring galaxy with its own 
central BH and both holes become gravitationally bound forming a binary. If the BHs 
are of comparable mass most of the angular momentum of the binary will reside in the 
orbital component of the angular momentum and, in general, this component will not 
be aligned with the spin of the AGN black hole. After the merger, the remnant resulting 
from the binary merger will have a final spin that will depend of the contribution of the 
binary's orbital angular momentum, the two pre-merger BH spins and the amount of 
angular momentum lost to gravitational waves. Figure 3 shows a diagram that reflects 
this interchange. This model was originally posed by Merritt and Ekers [33] and Zier 
and Biermann [34] and later confirmed by the full relativistic simulations of Campanelli 
et al. [35] . There it was shown that the change in direction between the spin of the larger 
hole and the spin of the final BH can exceed 90 degrees, consistent with observations 
of radio X-shaped patterns. They also exemplify the astrophysically relevant case of 
gradual semi-periodical shifts in the jet direction. Changes of this class are related to 
the precession of the BH spin during the inspiral and have also been observed in radio 
[36]. 

2.4. Mass and Spin of the Final Black Hole 

The merger of a BBH will form a final larger BH of mass Mf and spin Sf. Thus, the 
initial state is described by eight parameters: the mass ratio q, the spin components 
of the two pre-merger BHs and the binary's angular velocity u. Here u specifies the 
starting point of a possibly very long inspiral trajectory. After the merger the final 
BH is characterized by seven parameters, the final mass Mf, the spin Sf and the kick 
velocity k. Predicting the final mass and spin from the initial parameters is of great 
importance in many astrophysical merger scenarios. 
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Figure 3. Binary black liole merger model of the change in the direction of the spin 
observed in some supermassive BHs. Two galaxies merge forming a BBH (left) that 
in general will have BH spins and orbital angular momentum in arbitrary directions. 
The black hole of mass M is surrounded by an accretion disk that powers the jets that 
are originally aligned with the spin S. After the merger (right) the newly formed BH 
will have a new spin direction Sf that, if the original BHs are of comparable mass, will 
be mostly dictated by the pre-merger orbital angular momentum L. 



Several groups have developed formulas that attempt to predict the final spin of the 
merger. The analytic estimate of Buonanno et al. [37] can give the final spin magnitude 
to within few percent with larger deviations for spins close to anti- alignment. There 
are also formulas which have been derived by fitting to results from full numerical 
simulations [16, 38, 39]. The approach by Campanelli et al. [16] gives only the magnitude 
of the final spin. This approach was later improved by Lousto et al. [40] who give 
equations for all final spin components, the final mass and the kick velocity. However, 
the formulas in [40] still depend on the infall direction at merger, i.e. on a quantity that 
is not known without performing numerical simulations. The approach by Barausse 
and RezzoUa [39] gives only the final spin, but all fitting coefficients are explicitly given. 
The formulas by Tichy and Marronetti [38] predict the the final spin vector as well as 
the final mass, and also do not contain any unknown coefficients. Comparisons of these 
formulas [39, 41] seem to indicate that both the approaches in [38] and [39] give similar 
results and have comparable errors. However, these comparisons were done for low or 
moderate initial spins. Lovelace et al. [42] have performed two numerical simulations 
that test these different formulas for the case of two initial BHs with very high spin 
(dimensionless spin magnitudes above 0.94). In these two cases the formula of Tichy 
and Marronetti [38] predicts a final spin value that is closest to the true numerical 
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answer. However, more high spin cases with mass ratios different from one need to be 
studied before one can draw firm conclusions. 

These formulas are useful in astrophysical models where BBH mergers are 
important. E.g. one can study how successive mergers of BBH [38] influence the final 
spin. Using the formula in [38] for the final mass both O'Neill et al. [43] and Rossi 
et al. [44] have considered how the mass loss due to merger affects a circumbinary 
disk. O'Neill et al. argue that the reduction of luminosity caused by the retreat of the 
inner edge of the disk following mass loss could be detectable. The same scenario is 
investigated by Rossi et al. [44] who cast doubt on its detectability, unless the final BH 
receives a substantial kick directed close to the disk plane. 

3. Gravitational Wave Observatories and Numerical Relativity 

As the latest generation of gravitational wave detectors comes on-line, the problem of 
faithfully simulating the evolution of binary systems of compact objects has become 
increasingly important. The detectors (NSF's LIGO [45], VIRGO [46], TAMA 
[47], GEO600 [48]) use laser interferometry to measure the tiny strains associated 
with passing gravitational waves [49], offering much higher sensitivity than previous 
experiments. In addition, new detectors are under planning and construction stages 
such as the Einstein Telescope [50], Indigo [51] and DECIGO [52]. As we mentioned 
in Section 1, due to the very low signal-to- noise ratio LIGO/ VIRGO (L/V) data 
analysts employ a method known as matched filtering in order to boost the chances of 
detectability. This technique requires a priori knowledge of the shape of a given signal 
which is compared with observational data to assess the probability of its presence 
amidst the noise. LIGO/VIRGO scientists and collaborators have created a bank of 
templates based on post-Newtonian methods that cover the early stages of the binary 
life (inspiral) which have been successfully incorporated to the data analysis pipelines 
used in all the scientific runs corresponding to the initial stage of the detectors §. 

Now that the numerically relativity community has achieved the maturity level 
needed for BBH simulations, many groups in the field (including ours) have clustered 
around several collaborative projects that directly involve L/V data analysts. Binary 
black holes in quasi-circular orbits can be described by eight parameters (see Section 2.4). 
Additionally, in order to remove any orbital eccentricity [53, 54, 55, 56, 57] that might be 
present in the initial data, a careful determination of the initial tangential momentum 
ensues adding an extra parameter ||. Covering a nine- dimensional parameter space with 
month-long simulations is such a daunting task that is clear that the L/V needs can only 
be met (at least in the short and intermediate term) by a close collaboration between 
numerical and analytical relativists. This common understanding in the relativity 

§ Since October 2010, LIGO has been taken off-line to undergo the upgrades leading to the next 
generation observatory Advanced LIGO. 

II Gravitational radiation has a circularizing effect on the BBH orbits. Any eccentricity found in the 
BBH at the moment of formation is expected to have radiated away by the time of the merger. 
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community has given birth to the NRAR (Numerical Relativity - Analytical Relativity) 
collaboration [58]. In addition to this effort, numerical groups around the world are 
also converging into a tightly knit collaboration with L/V data analysts that aims at 
incorporating numerically generated gravitational waveforms into the observatories data 
pipelines: the NINJA collaboration [59]. 

3.1. NINJA and NRAR Collaborations 

The Numerical Injection Analysis (NINJA) project started more than two years ago and 
consists of the largest scientific collaboration between numerical relativists and L/V data 
analysts to date. NINJA's goal is the smooth incorporation of BBH numerical templates 
into the observatories detection infrastructure. During the first stage of the collaboration 
the waveforms provided by the numerical groups were embedded in colored Gaussian 
noise and injected into the data analysis pipelines. After these results were published 
[60, 61], plans for a follow-up project were drawn giving birth to NINJA-2. The goal 
of the second stage of the project is producing longer and more accurate gravitational 
waveforms by the numerical groups and executing more systematic tests, this time using 
real L/V data noise. Our team has been part of NINJA since its inception. NINJA- 
2 has set guidelines on the length and quality of the waveforms to be prepared by 
the numerical groups. One of them is that the simulations should cover at least ten 
waveform cycles (five BBH orbits) not including the short period at the beginning of 
the simulation when the signal is corrupted by the "junk" radiation found in the initial 
data sets. Other restrictions are that the error accumulated in the amplitude and phase 
of the waveforms I = 2 m = 2 mode during the inspiral should not exceed 5% and 0.5 
radians respectively. 

The Numerical Relativity - Analytical Relativity (NRAR) collaboration was kick- 
started at the end of 2009 and consists of a thrust by a large group (more than 50 
members) of numerical and analytical relativists. The goal is to coordinate recent 
developments in the production of BBH waveforms in both numerical and analytical 
methods. One of the goals of the collaboration is the generation of a small survey 
of high-accuracy relatively long (20 cycles or more) numerical waveforms to be used 
in the calibration of analytical template banks. These can then be used to generate 
efficiently and faithfully tens of thousands of templates for L/V matched-filtering search 
algorithms. Like in the case of NINJA-2, the NRAR collaboration has also set guidelines 
on the length and quality of the waveforms. NRAR requires waveforms with at least 
twenty waveform cycles and an accumulated phase error not to exceed 0.05 radians. In 
addition to these criteria, NRAR opted for minimizing the initial orbital eccentricity to 
make the signals more "astrophysically realistic" (see Section 2.4). 

3.2. Characteristics of a Typical Binary Black Hole Simulation 

In this section we describe some of the characteristics common to the most advanced 
simulations that are being produced by our group for the NINJA-2 and NRAR projects. 



Binary Black Holes: Recent Advances 



10 




Figure 4. The solid black lines mark the boundary of each refinement level. Shown 
are four concentric fixed and two moving levels around each BH. Each level doubles 
the spatial resolution. 
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Figure 5. Computing time (measured in wall clock days) vs. number of cores for 
a full-grid BBH simulation such as the one described in Section 3.2. Scaling was 
performed on NICS Cray XT-5 Kraken. 



The simulations are performed with the code BAM [62, 63, 64, 65, 25, 38] which evolves 
the gravitational fields using the BSSNOK formalism [66, 67, 68] in the variation known 
as the "moving punctures" method [11, 12]. BAM is based on a method of lines approach 
using sixth order finite differencing in space and explicit fourth order Runge-Kutta (RK) 
time stepping. The particulars of our numerical implementation can be found in [69, 65]. 
The numerical domain is represented by a hierarchy of nested Cartesian boxes as shown 
in a simplified manner in Fig. 4. It consists of L + 1 levels of refinement, indexed by 
I = 0, . . . ,L. A refinement level consists of one or two Cartesian boxes with a constant 
grid-spacing hi = ho/2'' on level /. Typical values are L + 1 = 12 for the number 
of refinement levels, with the levels through 5 each consisting of a single fixed box 
centered on the origin (the center of mass). On each of the finer levels (6 through L), 
we initially use two sets of moving boxes centered on each BH. When the BHs get close 
enough that two of these boxes start touching, they are replaced by a single box. The 
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Figure 6. Left: Gravitational wave strain with polarizations h^ (top) and hy. 
(bottom) times the distance to the observer D vs. time measured in units of total 
mass. The signal corresponds to an equal mass binary with spins of equal magnitude 
si/rn,1 — S2/m2 = 0.6 but arbitrary orientations. 



scaling of wall-clock time with number of processors is given in Fig. 5. 

Figure 2 shows the trajectory of the individual BHs in an equal mass binary. The 
initial BH spins are s1/mj = (0.5196, 0, 0.3000) and sl/m^ = (0, 0.5196, -0.3000) which 
have the magnitude si/mf = 82/1712 = 0.6. These parameters make the simulation fully 
three-dimensional: there are no underlying symmetries that can be exploited for the 
sake of computational efficiency. The system is evolved for about 10 orbits and the 
gravitational waveforms corresponding to the dominant quadrupole mode {I = m = 2) 
are shown in Fig. 6. Figure 7 compares the largest contribution for the multipoles with 
{1 = 2 and 4), highlighting the dominance of the quadrupole component. In both graphs 
the time is expressed in units of M, the total mass of the system. The merger produces 
a final BH with a total mass Mf = 0.952 (from an initial fiducial mass Mi = 1.0), with 
dimensionless spin of Sf/Mj = 0.704 and a recoil velocity of 398 km/ s. 



4. Summary 

We have reviewed some of the most impressive achievements produced by the numerical 
solution of the Einstein field equations corresponding to BBH. A long path has been 
traveled in the few years since the community finally came to master the intricacies of 
these simulations. A much longer road awaits ahead. Currently, most groups that used 
to work on vacuum problems are making inroads into simulations with matter such 
as binaries with neutrons stars and white dwarfs, core-collapse supernova and BBH 
embedded in circumbinary disks. To the challenges inherent to dealing with the rich 
microphysics required by the latter simulations, we should add numerical challenges such 
as preparing the code for the next generation of peta/exa-scale platforms. While these 
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Figure 7. Left: Gravitational wave strain with polarization h+ times the distance to 
the observer D vs. time measured in units of total mass. The curves shows the largest 
m modes for the moments I — 2, 4 (cases where / = m) of the simulations from Fig. 
6. It can be seen that most of the power is emitted in the quadrupole {I — 2) mode. 



problems are at least as challenging as the ones from the past, the results from these 
short years indicate that numerical relativity has a finally entered its golden age and, 
anticipating the detection of gravitational waves by the new generation of observatories, 
it couldn't happen at a better time. 
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